Name: Begum Babur UNI: bgb2123
Let’s load our packages first. NOTE: To install these packages in MacOS, you have to make sure GDAL and Gfortran are installed in your system. You can do this by brew install gdal and brew install gfortran via the Homebrew package manager.
# install.packages("rgdal")
# install.packages("fields")
# install.packages("dlm")
# install.packages("binhex")
library(jsonlite)
library(sp)
library(rgdal)
library(fields)
library(lubridate)
library(tidyverse)
library(ggplot2)
library(dlm)
Don’t forget to set your working directory, wherever your files are. NOTE: Below line will only work for RStudio users; change the argument to setwd() if you are not using RStudio.
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
Check whether you are in the current working directory. NOTE: Don’t forget to put the provided gps and test_gps folders, as well as the Weather.csv into the same directory this notebook is in.
list.files()
[1] "gps" "gps.zip" "Initial Analysis.R" "main.nb.html" "main.Rmd" "README.md"
[7] "test_gps" "test_gps.zip" "weather.csv"
Let’s create a function which will accepts a day’s .geojson and return a data frame that consists of: i) latitude, ii) longtitude, iii) time, iv) time differences, and v) speed.
create_df <- function(geojson_path) {
# Read in the .geojson file
stopifnot(endsWith(geojson_path, '.geojson'))
geojson = read_json(geojson_path, simplifyVector = T)
# Get the coordinates from the geojson
coordinates <- geojson[["features"]][["geometry"]][["coordinates"]]
# Split coordinates to latitude and longitude
latitude <- sapply(strsplit(as.character(coordinates)," "), "[[", 2)
latitude <- gsub("\\)", "", latitude)
latitude <- as.numeric(latitude)
longitude <- sapply(strsplit(as.character(coordinates)," "), "[[", 1)
longitude <- gsub("c\\(", "", longitude)
longitude <- gsub(",", "", longitude)
longitude <- as.numeric(longitude)
# Modify time vector and adjust for timezone
time <- geojson[["features"]][["properties"]][["time"]]
time <- as.POSIXct(strptime(time, tz = "MST", format = "%Y-%m-%dT%H:%M:%OS"))
time_long <- geojson[['features']]$properties$time_long
# Create data frame of coordinates and time values.
df <- data.frame(latitude, longitude, time, time_long)
# Include time differences between the recordings and the ones above.
df$time.diff.1 <-lag(lead(df$time, 1) - df$time)
# Add speed values.
df$speed <- geojson[["features"]][["properties"]][["speed"]]
return(df)
}
Let’s run the above function on the corresponding GEOJSON data for the each of the 11 days.
df1 <- create_df(geojson_path = "gps/20200818114606.geojson")
df2 <- create_df(geojson_path = "gps/20200819132607.geojson")
df3 <- create_df(geojson_path = "gps/20200820151044.geojson")
df4 <- create_df(geojson_path = "gps/20200821111447.geojson")
df5 <- create_df(geojson_path = "gps/20200824130857.geojson")
df6 <- create_df(geojson_path = "gps/20200825121346.geojson")
df7 <- create_df(geojson_path = "gps/20200826131614.geojson")
df8 <- create_df(geojson_path = "gps/20200827113234.geojson")
df9 <- create_df(geojson_path = "gps/20200828122627.geojson")
df10 <- create_df(geojson_path = "gps/20200828130816.geojson")
df11 <- create_df(geojson_path = "gps/20200831115147.geojson")
Let’s put the data frames together into one mega dataset which we can later use to fit some models!
# Combine all data frames in a list and change column names
df_list <- list(df1, df2, df3, df4, df5, df6, df7, df8, df9,df10, df11)
df_list <- lapply(df_list,
function(x) {
names(x) <- c("latitude", "longitude", "Time", "time_long", "time_difference", "speed")
return(x)
})
# Merge lists into one mega dataset
df <- do.call("rbind", df_list)
# Add combined coordinates to merged sets
df$coordinates <- paste(df$latitude, ",", df$longitude)
# Used this to get the most frequent coordinate, taking latitude and longitude values of each recording as a whole.
# Creating a mode function for coordinates in some sense.
all_coordinates <- table(df$coordinates)
Get the modes for the coordinates.
all_coordinates[all_coordinates == max(all_coordinates)]
46.8600804 , -113.9839519 46.8616699 , -113.9855465 46.8689186 , -113.9893428 46.8701026 , -113.9918634 46.8713271 , -113.9926369
3 3 3 3 3
46.8724778 , -113.9952979 46.8734358 , -113.9964836 46.8743452 , -113.9940895 46.8821394 , -113.998802 46.88427 , -113.9968601
3 3 3 3 3
df
Let’s get weather info into to add later to our dataset.
weather <- read.csv("weather.csv", stringsAsFactors = FALSE, header = TRUE)
weather <- weather[, 2:3]
# Format weather date and time
weather$Time <- as.POSIXct(lubridate::parse_date_time(weather$Date_Time, tz ="MST", orders='mdyhm'))
hms, hm and ms usage is deprecated, please use HMS, HM or MS instead. Deprecated in version '1.5.6'.
weather$Temp <- as.numeric(weather$Temp)
weather
Let’s convert the time formats of both datasets so that they match and can be left-joined later.
weather$Time <- format(weather$Time,"%D-%H")
weather
df$Time <- format(df$Time,"%D-%H")
df
Now, let’s merge the two data frames (df and weather) by their timestamps.
combined <- left_join(df, weather, by="Time")
combined
Let’s create a hex plot and show the frequencies of the coordinates (for both weeks).
plot_all <- ggplot(df, mapping = aes(latitude, longitude)) +
geom_hex() +
geom_point(x=46.8600804, y=-113.9839519, color = "red") + # Points created according to mode values
geom_point(x=46.8616699, y=-113.9855465, color = "red") +
geom_point(x=46.8689186, y=-113.9893428, color = "red") +
geom_point(x=46.8701026, y=-113.9918634, color = "red") +
geom_point(x=46.8713271, y=-113.9952979, color = "red") +
geom_point(x=46.8734358, y=-113.9964836, color = "red") +
geom_point(x=46.8743452, y=-113.9940895, color = "red") +
geom_point(x=46.8821394, y=-113.998802, color = "red") +
geom_point(x=46.88427, y=-113.9968601, color = "red")
# However, the modes of the repeated coordinates do not match the counts of the graph.
plot_all
Let’s now analyze the first week.
# Collect the data frames that belong to the 2nd week only
df_first_week_list <- list(df1, df2, df3, df4, df5)
df_first_week_list <- lapply(df_first_week_list,
function(x) {
names(x) <- c("latitude", "longitude", "time_MST","time_long", "time_difference", "speed")
return(x)
})
# Merge second week's list into one dataset
df_first_week <- do.call("rbind", df_first_week_list)
# Add coordinate to merged data frame
df_first_week$coordinates <- paste(df_first_week$latitude, ",", df_first_week$longitude)
all_coordinates_first_week <- table(df_first_week$coordinates)
Get the modes for the coordinates of the 2nd week only.
all_coordinates_first_week[all_coordinates_first_week == max(all_coordinates_first_week)]
46.8701026 , -113.9918634
3
Now let’s modify the time according to the day-hour format, and merge weather data frame with the data frame for the first week, just like what we previously did.
df_first_week$Time <- format(df_first_week$time_MST, "%D-%H")
combined_first_week <- left_join(df_first_week, weather, by="Time")
combined_first_week
Let’s create a hex plot for the first week only and show the frequencies of the coordinates. NOTE: The brightest hexagons do not correspond to the most repeated coordinates.
plot_first_week <- ggplot(df_first_week, mapping = aes(latitude, longitude)) +
geom_hex() +
geom_point(x=46.8701026, y=-113.9918634, color = "red")
plot_first_week
Let’s now analyze the second week.
# Collect the data frames that belong to the 2nd week only
df_second_week_list <- list(df6, df7, df8, df9, df10, df11)
df_second_week_list <- lapply(df_second_week_list,
function(x) {
names(x) <- c("latitude", "longitude", "time_MST","time_long", "time_difference", "speed")
return(x)
})
# Merge second week's list into one dataset
df_second_week <- do.call("rbind", df_second_week_list)
# Add coordinate to merged data frame
df_second_week$coordinates <- paste(df_second_week$latitude, ",", df_second_week$longitude)
all_coordinates_second_week <- table(df_second_week$coordinates)
Get the modes for the coordinates of the 2nd week only.
all_coordinates_second_week[all_coordinates_second_week == max(all_coordinates_second_week)]
46.8600804 , -113.9839519 46.8616699 , -113.9855465 46.8689186 , -113.9893428 46.8713271 , -113.9926369 46.8724778 , -113.9952979
3 3 3 3 3
46.8734358 , -113.9964836 46.8743452 , -113.9940895 46.8821394 , -113.998802 46.88427 , -113.9968601
3 3 3 3
Now let’s modify the time according to the day-hour format, and merge weather data frame with the data frame for the second week, just like what we previously did.
df_second_week$Time <- format(df_second_week$time_MST, "%D-%H")
combined_second_week <- left_join(df_second_week, weather, by="Time")
combined_second_week
Let’s create a hex plot for the second week only and show the frequencies of the coordinates. NOTE: The brightest hexagons do not correspond to the most repeated coordinates.
plot_second_week <- ggplot(df_second_week, mapping = aes(latitude, longitude)) +
geom_hex() +
geom_point(x=46.8600804, y=-113.9839519, color = "red") +
geom_point(x=46.8616699, y=-113.9855465, color = "red") +
geom_point(x=46.8689186, y=-113.9893428, color = "red") +
geom_point(x=46.8713271, y=-113.9926369, color = "red") +
geom_point(x=46.8724778, y=-113.9952979, color = "red") +
geom_point(x=46.8734358, y=-113.9964836, color = "red") +
geom_point(x=46.8743452, y=-113.9940895, color = "red") +
geom_point(x=46.8821394, y=-113.998802, color = "red") +
geom_point(x=46.88427, y=-113.9968601, color = "red")
plot_second_week
Let’s now go back to the data frame that contains both weeks, df. Remove the missing values and finding the median.
second_per_rec = median(df$time_difference, na.rm=TRUE)
second_per_rec
Time difference of 8 secs
The median time difference is 8 seconds. There are periods of inactivity, where the subject is stationary. These indicate times that we can’t bomb him.
print(paste("Number of seconds per record:", second_per_rec))
[1] "Number of seconds per record: 8"
print(paste("Number of records per second:", 1 / as.numeric(second_per_rec)))
[1] "Number of records per second: 0.125"
When the target is stationary for over 2 minutes, he’s considered safe. We will record this as the number of sessions of inactivity.
inactive <- length(df[df$time_difference > 120,])
inactive
[1] 7
In this dataset, he’s inactive 7 times.
Let’s apply the spatial transformation.
spat_df <- SpatialPointsDataFrame(coords = df[, c("longitude", "latitude")],
data = df[, c("longitude", "latitude", "Time")],
proj4string = CRS("+proj=longlat +datum=WGS84"))
spat_df
Let’s convert the longitude/latitude to UTM.
utm_df <- spTransform(spat_df, CRSobj = "+proj=utm +zone=12 +datum=WGS84 +units=m")
utm_coords <-coordinates(utm_df)
utm_df
utm_coords
longitude latitude
[1,] 271422.2 5196954
[2,] 271407.1 5196903
[3,] 271456.3 5196870
[4,] 271504.3 5196886
[5,] 271614.3 5196845
[6,] 271656.7 5196850
[7,] 271662.6 5196778
[8,] 271611.1 5196720
[9,] 271610.2 5196710
[10,] 271619.6 5196634
[11,] 271680.7 5196573
[12,] 271645.1 5196499
[13,] 271638.4 5196465
[14,] 271591.2 5196461
[15,] 271648.0 5196376
[16,] 271655.9 5196313
[17,] 271676.4 5196291
[18,] 271710.9 5196245
[19,] 271694.4 5196151
[20,] 271652.5 5196144
[21,] 271681.1 5196038
[22,] 271693.2 5196044
[23,] 271710.4 5196026
[24,] 271791.3 5195954
[25,] 271836.0 5195966
[26,] 271722.2 5195854
[27,] 271682.0 5195818
[28,] 271675.3 5195812
[29,] 271638.3 5195741
[30,] 271603.5 5195672
[31,] 271557.8 5195640
[32,] 271575.7 5195608
[33,] 271517.5 5195500
[34,] 271514.3 5195493
[35,] 271502.3 5195461
[36,] 271562.0 5195436
[37,] 271633.6 5195428
[38,] 271632.3 5195423
[39,] 271709.3 5195357
[40,] 271770.0 5195339
[41,] 271770.0 5195339
[42,] 271828.4 5195296
[43,] 271852.4 5195280
[44,] 271925.3 5195267
[45,] 271953.3 5195244
[46,] 271972.5 5195238
[47,] 272103.7 5195166
[48,] 272109.4 5195167
[49,] 272111.6 5195160
[50,] 272199.4 5195106
[51,] 272271.4 5195028
[52,] 272280.7 5195007
[53,] 272357.7 5195023
[54,] 272294.8 5194923
[55,] 272295.6 5194930
[56,] 272275.6 5194875
[57,] 272323.1 5194739
[58,] 272268.4 5194677
[59,] 272264.4 5194627
[60,] 272243.3 5194583
[61,] 272212.6 5194536
[62,] 272183.0 5194510
[63,] 272152.3 5194515
[64,] 272163.2 5194499
[65,] 272192.4 5194510
[66,] 272320.5 5194488
[67,] 272318.5 5194474
[68,] 272324.3 5194406
[69,] 272345.2 5194375
[70,] 272361.4 5194251
[71,] 272364.5 5194234
[72,] 272343.8 5194142
[73,] 272363.6 5194125
[74,] 272404.3 5194081
[75,] 272429.2 5194056
[76,] 272454.5 5194023
[77,] 272445.0 5194054
[78,] 272441.5 5194050
[79,] 272540.9 5193983
[80,] 272539.9 5193915
[81,] 272067.5 5195191
[82,] 272012.5 5195205
[83,] 271946.4 5195242
[84,] 271946.7 5195272
[85,] 271966.4 5195309
[86,] 272021.4 5195363
[87,] 271926.1 5195336
[88,] 271919.3 5195435
[89,] 271912.2 5195464
[90,] 271954.1 5195480
[91,] 271410.5 5196982
[92,] 271410.5 5196982
[93,] 271436.9 5196925
[94,] 271445.5 5196885
[95,] 271445.5 5196880
[96,] 271515.5 5196857
[97,] 271560.7 5196863
[98,] 271604.6 5196845
[99,] 271652.6 5196859
[100,] 271652.6 5196859
[101,] 271660.3 5196783
[102,] 271608.3 5196721
[103,] 271607.1 5196713
[104,] 271603.0 5196553
[105,] 271607.0 5196552
[106,] 271602.7 5196561
[107,] 271674.7 5196565
[108,] 271721.0 5196544
[109,] 271716.4 5196507
[110,] 271716.4 5196507
[111,] 271770.4 5196486
[112,] 271730.5 5196444
[113,] 271698.2 5196402
[114,] 271700.8 5196400
[115,] 271686.5 5196355
[116,] 271671.2 5196303
[117,] 271661.3 5196299
[118,] 271661.3 5196299
[119,] 271726.8 5196262
[120,] 271713.2 5196210
[121,] 271689.5 5196188
[122,] 271697.2 5196149
[123,] 271728.0 5196118
[124,] 271730.4 5196116
[125,] 271773.0 5196117
[126,] 271808.4 5196050
[127,] 271822.6 5196048
[128,] 271811.6 5196001
[129,] 271826.5 5195964
[130,] 271713.0 5195879
[131,] 271685.0 5195854
[132,] 271677.9 5195812
[133,] 271675.0 5195813
[134,] 271644.1 5195749
[135,] 271638.1 5195724
[136,] 271611.5 5195670
[137,] 271593.0 5195652
[138,] 271610.7 5195671
[139,] 271608.5 5195669
[140,] 271554.7 5195562
[141,] 271517.6 5195499
[142,] 271513.6 5195505
[143,] 271515.4 5195487
[144,] 271520.5 5195446
[145,] 271521.1 5195448
[146,] 271524.4 5195442
[147,] 271598.1 5195443
[148,] 271597.4 5195414
[149,] 271640.6 5195388
[150,] 271667.7 5195376
[151,] 271686.2 5195365
[152,] 271686.2 5195365
[153,] 271750.3 5195341
[154,] 271813.7 5195310
[155,] 271813.7 5195310
[156,] 271836.4 5195292
[157,] 271911.9 5195263
[158,] 271925.3 5195263
[159,] 271928.6 5195268
[160,] 271972.3 5195236
[161,] 272024.6 5195192
[162,] 272070.3 5195188
[163,] 272103.5 5195157
[164,] 272110.6 5195172
[165,] 272137.0 5195161
[166,] 272150.8 5195146
[167,] 272159.9 5195097
[168,] 272226.9 5195047
[169,] 272216.2 5195055
[170,] 272209.6 5195025
[171,] 272285.0 5195012
[172,] 272280.2 5195010
[173,] 272298.8 5195003
[174,] 272322.9 5194968
[175,] 272313.2 5194945
[176,] 272307.8 5194931
[177,] 272292.2 5194901
[178,] 272265.4 5194861
[179,] 272261.8 5194842
[180,] 272296.8 5194786
[181,] 272304.0 5194723
[182,] 272280.4 5194684
[183,] 272266.9 5194643
[184,] 272258.9 5194623
[185,] 272233.8 5194593
[186,] 272161.9 5194511
[187,] 272148.4 5194498
[188,] 272158.9 5194514
[189,] 272158.9 5194514
[190,] 272194.0 5194518
[191,] 272159.8 5194509
[192,] 272200.7 5194502
[193,] 272278.7 5194491
[194,] 272310.5 5194489
[195,] 272319.8 5194483
[196,] 272336.3 5194435
[197,] 272330.6 5194400
[198,] 272330.6 5194400
[199,] 272305.8 5194369
[200,] 272378.8 5194324
[201,] 272368.0 5194275
[202,] 272387.5 5194239
[203,] 272387.5 5194239
[204,] 272439.7 5194200
[205,] 272439.2 5194198
[206,] 272399.6 5194129
[207,] 272399.6 5194129
[208,] 272473.1 5194078
[209,] 272469.8 5194044
[210,] 272455.9 5194013
[211,] 272449.8 5194009
[212,] 272462.6 5193948
[213,] 272474.1 5193908
[214,] 272457.3 5193928
[215,] 272462.0 5193938
[216,] 272466.9 5193981
[217,] 272453.7 5194014
[218,] 272465.0 5194020
[219,] 272458.0 5194036
[220,] 272462.2 5194066
[221,] 272439.3 5194128
[222,] 272439.3 5194128
[223,] 272424.4 5194186
[224,] 272442.0 5194203
[225,] 272441.9 5194218
[226,] 272431.6 5194297
[227,] 272407.7 5194308
[228,] 272371.7 5194362
[229,] 272351.4 5194362
[230,] 272353.7 5194394
[231,] 272374.3 5194431
[232,] 272358.6 5194450
[233,] 272325.7 5194426
[234,] 272343.7 5194511
[235,] 272337.1 5194541
[236,] 272306.9 5194558
[237,] 272287.7 5194590
[238,] 272294.0 5194592
[239,] 272295.6 5194675
[240,] 272284.4 5194703
[241,] 272294.3 5194742
[242,] 272287.5 5194778
[243,] 272341.1 5194895
[244,] 272341.1 5194895
[245,] 272301.8 5194903
[246,] 272291.8 5194898
[247,] 272333.4 5194917
[248,] 272379.1 5194970
[249,] 272365.6 5194990
[250,] 272383.9 5195046
[251,] 272395.2 5195069
[252,] 272394.8 5195067
[253,] 272394.8 5195067
[254,] 272293.3 5195089
[255,] 272302.2 5195103
[256,] 272302.2 5195103
[257,] 272234.6 5195143
[258,] 272234.6 5195143
[259,] 272147.2 5195151
[260,] 272147.2 5195151
[261,] 272110.3 5195188
[262,] 272110.3 5195188
[263,] 272037.3 5195248
[264,] 271956.2 5195242
[265,] 271956.2 5195242
[266,] 271949.6 5195246
[267,] 271949.6 5195246
[268,] 271958.2 5195301
[269,] 271958.2 5195301
[270,] 272018.5 5195379
[271,] 272018.5 5195379
[272,] 271911.7 5195397
[273,] 271909.5 5195421
[274,] 271909.5 5195421
[275,] 271898.3 5195433
[276,] 271911.4 5195519
[277,] 271880.2 5195544
[278,] 271860.4 5195573
[279,] 271855.0 5195552
[280,] 271811.4 5195592
[281,] 271775.1 5195604
[282,] 271779.2 5195603
[283,] 271704.3 5195625
[284,] 271672.2 5195630
[285,] 271667.3 5195636
[286,] 271651.6 5195636
[287,] 271641.7 5195648
[288,] 271626.6 5195711
[289,] 271627.3 5195714
[290,] 271629.8 5195681
[291,] 271624.5 5195693
[292,] 271647.5 5195757
[293,] 271642.8 5195750
[294,] 271667.3 5195797
[295,] 271675.7 5195801
[296,] 271687.8 5195828
[297,] 271819.8 5195981
[298,] 271828.4 5195990
[299,] 271832.3 5195980
[300,] 271847.9 5196008
[301,] 271825.6 5196044
[302,] 271825.6 5196044
[303,] 271789.6 5196069
[304,] 271720.2 5196134
[305,] 271720.5 5196129
[306,] 271721.5 5196132
[307,] 271690.7 5196214
[308,] 271721.1 5196201
[309,] 271718.3 5196239
[310,] 271718.3 5196239
[311,] 271719.2 5196299
[312,] 271635.3 5196345
[313,] 271649.9 5196324
[314,] 271619.5 5196346
[315,] 271617.2 5196382
[316,] 271611.6 5196418
[317,] 271596.7 5196461
[318,] 271645.4 5196476
[319,] 271597.8 5196536
[320,] 271615.6 5196541
[321,] 271613.9 5196557
[322,] 271625.9 5196644
[323,] 271612.8 5196629
[324,] 271616.9 5196705
[325,] 271614.4 5196752
[326,] 271636.7 5196767
[327,] 271648.3 5196792
[328,] 271610.4 5196834
[329,] 271585.7 5196873
[330,] 271587.6 5196867
[331,] 271519.4 5196879
[332,] 271517.4 5196864
[333,] 271467.7 5196877
[334,] 271418.4 5196853
[335,] 271436.4 5196917
[336,] 271408.4 5196970
[337,] 271427.1 5196981
[338,] 271442.3 5197006
[339,] 271445.0 5197021
[340,] 271439.6 5196993
[341,] 271421.0 5196961
[342,] 271437.3 5196926
[343,] 271402.6 5196935
[344,] 271402.6 5196935
[345,] 271481.3 5196850
[346,] 271522.9 5196845
[347,] 271553.6 5196860
[348,] 271609.1 5196840
[349,] 271641.3 5196853
[350,] 271641.3 5196853
[351,] 271677.4 5196827
[352,] 271606.9 5196746
[353,] 271608.7 5196719
[354,] 271601.3 5196717
[355,] 271608.5 5196705
[356,] 271606.5 5196704
[357,] 271614.4 5196694
[358,] 271612.9 5196696
[359,] 271626.9 5196604
[360,] 271614.4 5196564
[361,] 271614.4 5196564
[362,] 271618.3 5196565
[363,] 271686.3 5196570
[364,] 271704.8 5196554
[365,] 271718.7 5196561
[366,] 271710.4 5196514
[367,] 271710.2 5196511
[368,] 271737.7 5196497
[369,] 271736.7 5196496
[370,] 271786.1 5196473
[371,] 271766.9 5196460
[372,] 271726.3 5196441
[373,] 271726.3 5196441
[374,] 271701.3 5196412
[375,] 271706.2 5196415
[376,] 271750.0 5196376
[377,] 271757.1 5196371
[378,] 271763.5 5196356
[379,] 271762.1 5196354
[380,] 271770.2 5196334
[381,] 271819.0 5196372
[382,] 271766.7 5196313
[383,] 271745.4 5196292
[384,] 271757.6 5196276
[385,] 271733.3 5196250
[386,] 271729.8 5196237
[387,] 271729.8 5196237
[388,] 271723.0 5196234
[389,] 271702.4 5196170
[390,] 271699.1 5196146
[391,] 271679.6 5196158
[392,] 271679.6 5196158
[393,] 271611.0 5196157
[394,] 271633.9 5196142
[395,] 271633.9 5196142
[396,] 271652.7 5196080
[397,] 271652.7 5196080
[398,] 271667.8 5196072
[399,] 271690.3 5196048
[400,] 271675.3 5196060
[401,] 271705.8 5196030
[402,] 271701.6 5196032
[403,] 271698.6 5196038
[404,] 271698.6 5196038
[405,] 271709.3 5196025
[406,] 271708.8 5196026
[407,] 271697.4 5196037
[408,] 271693.4 5196044
[409,] 271684.6 5196069
[410,] 271639.4 5195982
[411,] 271680.5 5196062
[412,] 271626.7 5195969
[413,] 271635.9 5195968
[414,] 271618.5 5195932
[415,] 271621.8 5195921
[416,] 271579.8 5195911
[417,] 271579.8 5195911
[418,] 271685.1 5195883
[419,] 271685.2 5195882
[420,] 271687.8 5195878
[421,] 271678.5 5195826
[422,] 271675.5 5195819
[423,] 271675.9 5195813
[424,] 271675.9 5195813
[425,] 271650.6 5195780
[426,] 271651.8 5195784
[427,] 271651.9 5195775
[428,] 271641.4 5195744
[429,] 271626.7 5195707
[430,] 271626.4 5195707
[431,] 271629.2 5195682
[432,] 271662.2 5195642
[433,] 271662.2 5195642
[434,] 271637.2 5195657
[435,] 271637.2 5195657
[436,] 271644.7 5195642
[437,] 271644.7 5195642
[438,] 271645.2 5195641
[439,] 271645.2 5195641
[440,] 271658.5 5195633
[441,] 271658.5 5195633
[442,] 271682.8 5195619
[443,] 271682.8 5195619
[444,] 271727.5 5195581
[445,] 271729.6 5195552
[446,] 271733.4 5195538
[447,] 271744.6 5195520
[448,] 271757.9 5195491
[449,] 271765.7 5195489
[450,] 271779.9 5195443
[451,] 271807.8 5195430
[452,] 271833.0 5195429
[453,] 271826.9 5195419
[454,] 271830.5 5195412
[455,] 271811.7 5195378
[456,] 271787.8 5195369
[457,] 271780.8 5195333
[458,] 271761.5 5195330
[459,] 271770.8 5195321
[460,] 271831.1 5195298
[461,] 271825.7 5195303
[462,] 271825.7 5195303
[463,] 271814.5 5195296
[464,] 271803.1 5195252
[465,] 271797.6 5195247
[466,] 271779.2 5195179
[467,] 271771.6 5195205
[468,] 271771.6 5195205
[469,] 271767.1 5195170
[470,] 271742.8 5195151
[471,] 271740.0 5195139
[472,] 271745.5 5195102
[473,] 271724.7 5195081
[474,] 271734.3 5195075
[475,] 271739.3 5195070
[476,] 271730.7 5195077
[477,] 271739.2 5195073
[478,] 271739.2 5195073
[479,] 271775.7 5195065
[480,] 271775.7 5195065
[481,] 271832.4 5195047
[482,] 271789.5 5194994
[483,] 271846.2 5195050
[484,] 271860.9 5195040
[485,] 271897.0 5195043
[486,] 271921.1 5195006
[487,] 271946.2 5195003
[488,] 271948.6 5194996
[489,] 271948.6 5194996
[490,] 271984.1 5194997
[491,] 272029.4 5194994
[492,] 272020.1 5194988
[493,] 272020.1 5194988
[494,] 272064.0 5194966
[495,] 272064.0 5194966
[496,] 272068.4 5194963
[497,] 272129.2 5194981
[498,] 272105.7 5194961
[499,] 272148.4 5194956
[500,] 272177.2 5194950
[ reached getOption("max.print") -- omitted 5357 rows ]
Let’s find the dates when he is stationary for more than 2 minutes.
which(df$time_difference > 120)
[1] 81 214 565 1022 1723 2404 3178 3535 3984 4661 5353 5519 5529
for (i in which(df$time_difference > 120)) {
print(df[i,]$Time)
}
[1] "08/18/20-22"
[1] "08/19/20-23"
[1] "08/21/20-00"
[1] "08/21/20-23"
[1] "08/24/20-22"
[1] "08/25/20-21"
[1] "08/26/20-23"
[1] "08/27/20-04"
[1] "08/27/20-22"
[1] "08/28/20-23"
[1] "08/31/20-19"
[1] "08/31/20-20"
[1] "08/31/20-22"
Aside from August 31st and 26th, he’s stationary at the night when he’s back home. Let’s drop these two dates.
df[5353,]$Time
[1] "08/31/20-19"
df[5529,]$Time
[1] "08/31/20-22"
GPS data reveals that on August 31st he’s at the University of Montana at that time. Removing these times when he’s at the university because he can’t be bombed when stationary.
nrow(df)
[1] 5857
df <- df[!(grepl('08/31/20', df$Time, fixed=TRUE) | grepl('08/26/20', df$Time, fixed=TRUE)),]
nrow(df)
[1] 4490
nrow(spat_df)
[1] 5857
spat_df <- spat_df[!(grepl('08/31/20', spat_df$Time, fixed=TRUE) | grepl('08/26/20', spat_df$Time, fixed=TRUE)),]
nrow(spat_df)
[1] 4490
Now let’s also spatially transform the week 2 dataset.
spat_df_second_week <-SpatialPointsDataFrame(coords = df_second_week[,c("longitude", "latitude")],
data = df_second_week[,c("longitude", "latitude", "Time")],
proj4string=CRS("+proj=longlat +datum=WGS84 +units=m"))
spat_df_second_week
Let’s converts the longitude/latitude to UTM
utm_df_second_week <-spTransform(spat_df_second_week, CRSobj = "+proj=utm +zone=12 +datum=WGS84 +units=m")
utm_coords_second_week <-coordinates(utm_df_second_week)
utm_df_second_week
utm_coords_second_week
longitude latitude
[1,] 271440.9 5197033
[2,] 271434.7 5197021
[3,] 271458.1 5197010
[4,] 271453.5 5196996
[5,] 271405.6 5196973
[6,] 271418.9 5196975
[7,] 271424.9 5196959
[8,] 271432.2 5196939
[9,] 271398.6 5196946
[10,] 271435.4 5196900
[11,] 271398.6 5196946
[12,] 271429.8 5196882
[13,] 271436.0 5196881
[14,] 271444.8 5196871
[15,] 271454.2 5196868
[16,] 271465.7 5196867
[17,] 271476.1 5196865
[18,] 271486.1 5196864
[19,] 271497.2 5196862
[20,] 271508.6 5196861
[21,] 271519.7 5196861
[22,] 271532.3 5196860
[23,] 271543.5 5196861
[24,] 271554.9 5196860
[25,] 271567.2 5196860
[26,] 271661.2 5196868
[27,] 271583.9 5196859
[28,] 271647.7 5196849
[29,] 271597.0 5196851
[30,] 271653.3 5196833
[31,] 271608.7 5196828
[32,] 271675.4 5196773
[33,] 271610.9 5196808
[34,] 271610.7 5196801
[35,] 271609.4 5196795
[36,] 271609.0 5196789
[37,] 271609.8 5196781
[38,] 271611.4 5196772
[39,] 271612.6 5196762
[40,] 271600.6 5196720
[41,] 271623.0 5196737
[42,] 271605.0 5196718
[43,] 271634.7 5196717
[44,] 271605.6 5196710
[45,] 271633.6 5196700
[46,] 271602.4 5196708
[47,] 271630.3 5196682
[48,] 271629.6 5196674
[49,] 271626.4 5196665
[50,] 271613.0 5196629
[51,] 271632.2 5196641
[52,] 271614.4 5196629
[53,] 271631.0 5196623
[54,] 271628.0 5196614
[55,] 271624.9 5196605
[56,] 271621.6 5196596
[57,] 271620.5 5196585
[58,] 271607.5 5196546
[59,] 271620.8 5196567
[60,] 271602.8 5196551
[61,] 271625.5 5196543
[62,] 271628.5 5196552
[63,] 271626.3 5196521
[64,] 271625.9 5196514
[65,] 271626.4 5196502
[66,] 271630.6 5196490
[67,] 271638.1 5196483
[68,] 271671.3 5196455
[69,] 271642.0 5196466
[70,] 271655.9 5196466
[71,] 271647.5 5196459
[72,] 271665.1 5196456
[73,] 271655.8 5196447
[74,] 271677.8 5196435
[75,] 271666.3 5196433
[76,] 271681.8 5196426
[77,] 271687.2 5196420
[78,] 271681.3 5196428
[79,] 271702.4 5196408
[80,] 271681.3 5196428
[81,] 271703.1 5196391
[82,] 271696.8 5196403
[83,] 271680.6 5196375
[84,] 271696.8 5196403
[85,] 271675.1 5196353
[86,] 271645.9 5196322
[87,] 271650.5 5196324
[88,] 271645.9 5196322
[89,] 271664.5 5196300
[90,] 271671.2 5196304
[91,] 271691.0 5196289
[92,] 271695.9 5196278
[93,] 271694.9 5196267
[94,] 271729.2 5196269
[95,] 271735.8 5196249
[96,] 271729.2 5196269
[97,] 271767.7 5196235
[98,] 271760.1 5196232
[99,] 271797.2 5196232
[100,] 271784.8 5196214
[101,] 271795.9 5196209
[102,] 271784.8 5196214
[103,] 271805.6 5196200
[104,] 271803.3 5196200
[105,] 271816.9 5196171
[106,] 271804.2 5196137
[107,] 271806.1 5196141
[108,] 271804.2 5196137
[109,] 271774.2 5196125
[110,] 271777.7 5196095
[111,] 271753.4 5196121
[112,] 271782.2 5196119
[113,] 271722.0 5196070
[114,] 271711.3 5196056
[115,] 271722.0 5196070
[116,] 271716.4 5196041
[117,] 271701.4 5196036
[118,] 271728.2 5196016
[119,] 271701.4 5196036
[120,] 271748.2 5196004
[121,] 271744.2 5196014
[122,] 271762.9 5195990
[123,] 271781.8 5195982
[124,] 271786.9 5195979
[125,] 271726.8 5195831
[126,] 271709.3 5195831
[127,] 271700.7 5195837
[128,] 271689.2 5195835
[129,] 271680.7 5195826
[130,] 271676.2 5195806
[131,] 271678.9 5195802
[132,] 271672.9 5195810
[133,] 271663.4 5195772
[134,] 271672.9 5195810
[135,] 271657.6 5195773
[136,] 271653.2 5195768
[137,] 271639.4 5195751
[138,] 271634.0 5195728
[139,] 271629.8 5195717
[140,] 271623.2 5195704
[141,] 271618.5 5195685
[142,] 271623.2 5195704
[143,] 271610.1 5195677
[144,] 271607.5 5195674
[145,] 271606.2 5195667
[146,] 271603.5 5195670
[147,] 271599.5 5195671
[148,] 271598.8 5195671
[149,] 271601.3 5195669
[150,] 271611.5 5195663
[151,] 271628.4 5195654
[152,] 271647.5 5195646
[153,] 271655.2 5195645
[154,] 271664.3 5195637
[155,] 271680.0 5195626
[156,] 271668.8 5195636
[157,] 271727.7 5195643
[158,] 271739.3 5195633
[159,] 271743.4 5195619
[160,] 271767.5 5195629
[161,] 271749.0 5195584
[162,] 271769.9 5195587
[163,] 271783.1 5195583
[164,] 271796.6 5195582
[165,] 271813.2 5195592
[166,] 271796.6 5195582
[167,] 271818.3 5195560
[168,] 271815.5 5195515
[169,] 271810.0 5195527
[170,] 271811.2 5195521
[171,] 271793.2 5195510
[172,] 271789.4 5195509
[173,] 271787.4 5195485
[174,] 271789.4 5195509
[175,] 271793.0 5195467
[176,] 271781.1 5195450
[177,] 271803.6 5195426
[178,] 271781.1 5195450
[179,] 271812.5 5195420
[180,] 271814.3 5195409
[181,] 271810.2 5195392
[182,] 271810.4 5195375
[183,] 271819.5 5195359
[184,] 271807.9 5195362
[185,] 271809.3 5195341
[186,] 271800.7 5195321
[187,] 271810.1 5195319
[188,] 271837.0 5195304
[189,] 271831.5 5195302
[190,] 271837.0 5195304
[191,] 271828.3 5195303
[192,] 271838.6 5195291
[193,] 271827.5 5195307
[194,] 271852.1 5195286
[195,] 271827.8 5195304
[196,] 271852.1 5195286
[197,] 271862.0 5195304
[198,] 271869.1 5195286
[199,] 271891.5 5195296
[200,] 271901.5 5195276
[201,] 271889.2 5195257
[202,] 271901.5 5195276
[203,] 271914.8 5195257
[204,] 271915.0 5195258
[205,] 271891.7 5195274
[206,] 271925.7 5195238
[207,] 271928.9 5195200
[208,] 271925.7 5195238
[209,] 271918.9 5195182
[210,] 271912.2 5195172
[211,] 271906.0 5195159
[212,] 271906.6 5195172
[213,] 271927.3 5195145
[214,] 271941.4 5195144
[215,] 271953.3 5195125
[216,] 271941.4 5195144
[217,] 271976.0 5195117
[218,] 271984.5 5195125
[219,] 272002.5 5195111
[220,] 272011.6 5195108
[221,] 272010.3 5195102
[222,] 271992.6 5195077
[223,] 272004.9 5195072
[224,] 271992.6 5195077
[225,] 272039.8 5195074
[226,] 272000.0 5195038
[227,] 272011.7 5195028
[228,] 271996.8 5195038
[229,] 272035.5 5194994
[230,] 272014.3 5194988
[231,] 272013.5 5194994
[232,] 272014.3 5194988
[233,] 272044.6 5194980
[234,] 272057.7 5194981
[235,] 272076.7 5194981
[236,] 272066.5 5194977
[237,] 272101.7 5194992
[238,] 272066.5 5194977
[239,] 272135.4 5194965
[240,] 272146.1 5194968
[241,] 272147.6 5194966
[242,] 272151.3 5194960
[243,] 272169.1 5194950
[244,] 272193.5 5194988
[245,] 272195.8 5194951
[246,] 272203.0 5194974
[247,] 272207.3 5194946
[248,] 272217.4 5194948
[249,] 272237.8 5194940
[250,] 272217.4 5194948
[251,] 272261.0 5194935
[252,] 272248.8 5194931
[253,] 272272.3 5194922
[254,] 272291.6 5194918
[255,] 272307.5 5194914
[256,] 272325.3 5194901
[257,] 272299.5 5194913
[258,] 272343.0 5194905
[259,] 272298.6 5194877
[260,] 272285.6 5194885
[261,] 272303.0 5194843
[262,] 272285.6 5194885
[263,] 272292.2 5194824
[264,] 272271.8 5194831
[265,] 272290.5 5194798
[266,] 272270.3 5194831
[267,] 272280.1 5194767
[268,] 272277.9 5194758
[269,] 272280.1 5194755
[270,] 272283.3 5194744
[271,] 272286.9 5194725
[272,] 272290.0 5194709
[273,] 272286.2 5194693
[274,] 272275.9 5194683
[275,] 272279.7 5194671
[276,] 272279.1 5194655
[277,] 272275.1 5194638
[278,] 272273.0 5194628
[279,] 272265.2 5194610
[280,] 272260.2 5194595
[281,] 272258.9 5194580
[282,] 272256.0 5194569
[283,] 272246.2 5194563
[284,] 272239.4 5194563
[285,] 272174.4 5194512
[286,] 272222.9 5194547
[287,] 272167.5 5194502
[288,] 272186.8 5194549
[289,] 272167.5 5194502
[290,] 272184.4 5194506
[291,] 272156.6 5194507
[292,] 272162.1 5194494
[293,] 272151.4 5194507
[294,] 272164.2 5194490
[295,] 272159.0 5194507
[296,] 272158.3 5194495
[297,] 272160.5 5194491
[298,] 272161.3 5194489
[299,] 272163.4 5194488
[300,] 272168.1 5194485
[301,] 272178.9 5194501
[302,] 272175.7 5194455
[303,] 272161.1 5194491
[304,] 272179.9 5194424
[305,] 272179.7 5194405
[306,] 272183.4 5194392
[307,] 272194.5 5194394
[308,] 272209.7 5194397
[309,] 272225.7 5194400
[310,] 272243.2 5194394
[311,] 272225.3 5194409
[312,] 272251.8 5194388
[313,] 272289.8 5194383
[314,] 272271.2 5194385
[315,] 272297.5 5194360
[316,] 272307.3 5194380
[317,] 272297.5 5194360
[318,] 272339.8 5194369
[319,] 272323.8 5194401
[320,] 272350.6 5194365
[321,] 272343.1 5194363
[322,] 272360.0 5194345
[323,] 272343.1 5194363
[324,] 272361.6 5194329
[325,] 272361.4 5194312
[326,] 272362.8 5194295
[327,] 272360.7 5194276
[328,] 272370.6 5194265
[329,] 272374.4 5194265
[330,] 272385.5 5194252
[331,] 272390.9 5194244
[332,] 272439.5 5194239
[333,] 272419.6 5194217
[334,] 272439.5 5194239
[335,] 272426.2 5194213
[336,] 272439.6 5194223
[337,] 272433.1 5194201
[338,] 272428.7 5194180
[339,] 272428.2 5194159
[340,] 272429.1 5194145
[341,] 272429.9 5194121
[342,] 272424.0 5194138
[343,] 272439.7 5194098
[344,] 272411.8 5194124
[345,] 272475.6 5194091
[346,] 272411.8 5194124
[347,] 272473.6 5194075
[348,] 272469.1 5194063
[349,] 272463.4 5194043
[350,] 272465.4 5194036
[351,] 272461.3 5194038
[352,] 272469.3 5194025
[353,] 272470.1 5194006
[354,] 272471.9 5193989
[355,] 272471.6 5193971
[356,] 272471.7 5193959
[357,] 272462.8 5193964
[358,] 272460.8 5193942
[359,] 272462.5 5193922
[360,] 272455.6 5194007
[361,] 272455.6 5194007
[362,] 272462.8 5193923
[363,] 272463.7 5193931
[364,] 272456.5 5194012
[365,] 272457.4 5193957
[366,] 272459.9 5193979
[367,] 272454.9 5193987
[368,] 272455.8 5193998
[369,] 272461.5 5194013
[370,] 272450.6 5194019
[371,] 272462.2 5194035
[372,] 272463.4 5194018
[373,] 272468.0 5194025
[374,] 272466.8 5194029
[375,] 272476.3 5194061
[376,] 272466.8 5194029
[377,] 272461.2 5194070
[378,] 272453.1 5194089
[379,] 272449.9 5194110
[380,] 272433.5 5194130
[381,] 272432.2 5194126
[382,] 272433.5 5194130
[383,] 272435.6 5194149
[384,] 272416.2 5194130
[385,] 272431.8 5194169
[386,] 272430.0 5194178
[387,] 272427.8 5194195
[388,] 272439.6 5194202
[389,] 272412.6 5194224
[390,] 272439.6 5194202
[391,] 272412.8 5194206
[392,] 272440.0 5194232
[393,] 272415.1 5194240
[394,] 272403.0 5194248
[395,] 272392.4 5194259
[396,] 272382.8 5194270
[397,] 272376.3 5194282
[398,] 272367.1 5194288
[399,] 272356.5 5194292
[400,] 272358.4 5194292
[401,] 272364.0 5194314
[402,] 272365.0 5194322
[403,] 272367.5 5194326
[404,] 272372.6 5194341
[405,] 272325.6 5194395
[406,] 272353.3 5194371
[407,] 272329.6 5194395
[408,] 272357.4 5194389
[409,] 272320.5 5194421
[410,] 272360.6 5194417
[411,] 272320.5 5194421
[412,] 272350.8 5194437
[413,] 272352.3 5194454
[414,] 272357.7 5194463
[415,] 272363.0 5194474
[416,] 272363.2 5194494
[417,] 272360.6 5194504
[418,] 272350.1 5194519
[419,] 272338.3 5194528
[420,] 272333.1 5194532
[421,] 272325.1 5194546
[422,] 272317.5 5194553
[423,] 272305.4 5194558
[424,] 272289.5 5194565
[425,] 272298.0 5194570
[426,] 272292.3 5194577
[427,] 272285.2 5194591
[428,] 272281.7 5194615
[429,] 272284.7 5194632
[430,] 272286.1 5194642
[431,] 272286.2 5194657
[432,] 272290.1 5194665
[433,] 272291.4 5194677
[434,] 272293.0 5194696
[435,] 272294.0 5194708
[436,] 272298.4 5194716
[437,] 272296.4 5194725
[438,] 272297.2 5194738
[439,] 272307.6 5194744
[440,] 272339.8 5194896
[441,] 272316.8 5194764
[442,] 272339.8 5194896
[443,] 272310.6 5194798
[444,] 272341.4 5194903
[445,] 272309.8 5194828
[446,] 272314.9 5194842
[447,] 272321.4 5194858
[448,] 272305.5 5194924
[449,] 272326.8 5194875
[450,] 272305.5 5194924
[451,] 272332.1 5194890
[452,] 272329.3 5194923
[453,] 272346.2 5194909
[454,] 272329.3 5194923
[455,] 272338.7 5194931
[456,] 272341.4 5194905
[457,] 272349.5 5194955
[458,] 272329.5 5194940
[459,] 272353.0 5194973
[460,] 272325.0 5194958
[461,] 272368.5 5195012
[462,] 272369.2 5195000
[463,] 272373.8 5195010
[464,] 272377.9 5195016
[465,] 272359.3 5195039
[466,] 272408.6 5195083
[467,] 272324.2 5195057
[468,] 272315.1 5195061
[469,] 272303.5 5195068
[470,] 272303.4 5195108
[471,] 272286.1 5195081
[472,] 272303.4 5195108
[473,] 272277.5 5195080
[474,] 272261.4 5195087
[475,] 272249.0 5195093
[476,] 272234.6 5195102
[477,] 272211.3 5195131
[478,] 272221.5 5195108
[479,] 272225.1 5195142
[480,] 272204.8 5195121
[481,] 272225.1 5195142
[482,] 272185.8 5195132
[483,] 272201.3 5195150
[484,] 272159.9 5195131
[485,] 272153.6 5195140
[486,] 272139.1 5195150
[487,] 272127.9 5195156
[488,] 272113.7 5195156
[489,] 272120.9 5195170
[490,] 272102.1 5195158
[491,] 272100.8 5195157
[492,] 272090.8 5195167
[493,] 272099.3 5195192
[494,] 272062.2 5195179
[495,] 272099.3 5195192
[496,] 272051.4 5195188
[497,] 272037.5 5195191
[498,] 272024.8 5195196
[499,] 272024.0 5195205
[500,] 272011.9 5195210
[ reached getOption("max.print") -- omitted 3313 rows ]
plot(coordinates(spat_df), col=factor(spat_df$Time))
Let’s manually create inbound and outbound trips according to your coordinates. 0 = inbound (coming home from work), 1 = outbound (leaving home for work)
spat_df$session <- rep(NA,nrow(spat_df@data))
spat_df$session[1:80] <- 0
spat_df$session[81:90] <- 1
spat_df$session[91:213] <- 0
spat_df$session[214:339] <- 1
spat_df$session[340:564] <- 0
spat_df$session[565:693] <- 1
spat_df$session[694:1021] <- 0
spat_df$session[1022:1360] <- 1
spat_df$session[1361:1722] <- 0
spat_df$session[1723:2044] <- 1
spat_df$session[2045:2403] <- 0
spat_df$session[2404:2755] <- 1
spat_df$session[2756:3117] <- 0
spat_df$session[3118:3533] <- 1
spat_df$session[3534:3896] <- 0
spat_df$session[3897:4245] <- 1
# spat_df$session[4246:4512] <- 0
# spat_df$session[4513:4573] <- 1
# spat_df$session[4574:4903] <- 0
# spat_df$session[4904:5264] <- 1
# spat_df$session[5265:5352] <- 0
outbound <- spat_df[spat_df[!is.na(spat_df$session==0),]$session==0,]
inbound <- spat_df[spat_df[!is.na(spat_df$session==1),]$session==1,]
outbound
inbound
Let’s plot the inbound and outbound plots separately.
plot(outbound$longitude, outbound$latitude)
points(inbound$longitude, inbound$latitude, col = "red")
gps_variance <- 20^2
v_mat <-diag(c(gps_variance, gps_variance))
f_mat <-matrix(c(1,0,0,0, 0,1,0,0), nrow=2, byrow = TRUE)
dt <- 8
g_mat <-matrix(c(1, 0, dt, 0,
0, 1, 0, dt,
0, 0, 1, 0,
0, 0, 0, 1), byrow=TRUE, ncol=4)
avg_walk_speed_m_per_sec <- median(df$speed, na.rm=TRUE)
dlm_spec <-dlm(FF= f_mat,
GG= g_mat,
V = v_mat,
W =diag(c(5, 5, 1, 1)^2),
m0 =matrix(c(utm_coords[1, ],rep(avg_walk_speed_m_per_sec/8, 2)),ncol=1),
# A vector by R defaults is a k by 1 matrix.
C0 =diag(rep(10^2, 4)))
dlm_filter_mod <- dlmFilter(utm_coords, dlm_spec)
dlm_smooth_mod <- dlmSmooth(dlm_filter_mod)
# Plot Kalman Filter, Smoother and Raw data into one plot to see how they're doing.
# See the first 100 values for a clearer picture.
par(mfrow=c(1,1))
plot(cbind(coordinates(spat_df)[1:200, ],
dlm_filter_mod$m[1:200, 1:2],
dlm_smooth_mod$s[1:200, 1:2]),ype='p', col =c("black", "red", "blue"), xlab="UTM X", ylab="UTM Y")
legend("topright", col = c("black", "red", "blue"), pch = 1,legend = c("raw", "kalman filter", "kalman smoother"))
axis(1, xaxp=c(46.85,46.9,10))
axis(2, yaxp=c(-114.07,-113.98,10))
Now let’s compare to other models. Using GLM, let’s predict latitude given longitude and time.
Let’s first look at the unique dates to make a decision on how we should form train and test splits.
unique(df_first_week$Time)
[1] "08/18/20-17" "08/18/20-18" "08/18/20-22" "08/19/20-19" "08/19/20-20" "08/19/20-23" "08/20/20-00" "08/20/20-21" "08/20/20-22" "08/21/20-00"
[11] "08/21/20-17" "08/21/20-18" "08/21/20-23" "08/22/20-00" "08/24/20-19" "08/24/20-22" "08/24/20-23"
unique(df_second_week$Time)
[1] "08/25/20-18" "08/25/20-21" "08/25/20-22" "08/26/20-19" "08/26/20-20" "08/26/20-23" "08/27/20-00" "08/27/20-04" "08/27/20-17" "08/27/20-18"
[11] "08/27/20-22" "08/27/20-23" "08/28/20-18" "08/28/20-19" "08/28/20-23" "08/31/20-17" "08/31/20-18" "08/31/20-19" "08/31/20-20" "08/31/20-22"
glm_latitude <- glm(data=df_first_week, latitude~longitude+time_long+speed+time_difference)
glm_predictions <- predict(glm_latitude, newdata=df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),]$latitude - glm_predictions
plot(glm_residuals)
glm_longitude <- glm(data=df_first_week, longitude~latitude+time_long+speed+time_difference)
glm_predictions <- predict(glm_longitude, newdata=df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),]$longitude - glm_predictions
plot(glm_residuals)
glm_time_long <- glm(data=df_first_week, time_long~longitude+latitude+speed+time_difference)
glm_predictions <- predict(glm_time_long, newdata=df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/25/20', df_second_week$Time, fixed=TRUE),]$time_long - glm_predictions
plot(glm_residuals)
glm_latitude <- glm(data=df_second_week[!grepl('08/31/20', df_second_week$Time, fixed=TRUE),], latitude~longitude+time_long+speed+time_difference)
glm_predictions <- predict(glm_latitude, newdata=df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),]$latitude - glm_predictions
plot(glm_residuals)
glm_longitude <- glm(data=df_second_week[!grepl('08/31/20', df_second_week$Time, fixed=TRUE),], longitude~latitude+time_long+speed+time_difference)
glm_predictions <- predict(glm_longitude, newdata=df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),]$longitude - glm_predictions
plot(glm_residuals)
glm_time_long <- glm(data=df_second_week[!grepl('08/31/20', df_second_week$Time, fixed=TRUE),], time_long~longitude+latitude+speed+time_difference)
glm_predictions <- predict(glm_time_long, newdata=df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),])
glm_residuals <- df_second_week[grepl('08/31/20', df_second_week$Time, fixed=TRUE),]$time_long - glm_predictions
plot(glm_residuals)
Now let’s code up our final algorithm with a single function! We will train/fit our models on the second week. k argument indicates the number of bombs you want to predict, based on the absolute sum error of residuals.
assasination_algorithm <- function(test_df, k=2) {
# Latitude prediction
glm_latitude <- glm(data=df_second_week, latitude~longitude+time_long+speed+time_difference)
glm_latitude_predictions <- predict(glm_latitude, newdata=test_df)
glm_latitude_residuals <- test_df$latitude - glm_latitude_predictions
# Longitude prediction
glm_longitude <- glm(data=df_second_week, longitude~latitude+time_long+speed+time_difference)
glm_longitude_predictions <- predict(glm_longitude, newdata=test_df)
glm_longitude_residuals <- test_df$longitude - glm_longitude_predictions
# Time prediction
glm_time_long <- glm(data=df_second_week, time_long~latitude+longitude+speed+time_difference)
glm_time_long_predictions <- predict(glm_time_long, newdata=test_df)
glm_time_long_residuals <- test_df$time_long - glm_time_long_predictions
# Measure errors (based on latitude and longitude only)
errors <- bind_cols(latitude=as.vector(glm_latitude_residuals), longitude=as.vector(glm_longitude_residuals))
absolute_sum_errors <- abs(errors$latitude) + abs(errors$longitude)
# Get bomb predictions (latitude, longitude, and time) based on top-k minimum errors
for (i in 1:k) {
print(paste('Bomb Prediction ', i))
print(paste('Latitude: ', glm_latitude_predictions[which(absolute_sum_errors == sort(absolute_sum_errors)[i])]))
print(paste('Longitude: ', glm_longitude_predictions[which(absolute_sum_errors == sort(absolute_sum_errors)[i])]))
print(paste('Time: ', glm_time_long_predictions[which(absolute_sum_errors == sort(absolute_sum_errors)[i])]))
print('-----------------------------------------------------------------------------------------------------------')
}
}
Now let’s read in the test data and format it accordingly, i.e. in the same way we have handled the training data.
test_df1 <- create_df(geojson_path = "test_gps/20200901112100.geojson")
test_df2 <- create_df(geojson_path = "test_gps/20200902125611.geojson")
test_df3 <- create_df(geojson_path = "test_gps/20200903110618.geojson")
test_df4 <- create_df(geojson_path = "test_gps/20200908081420.geojson")
test_df5 <- create_df(geojson_path = "test_gps/20200910070926.geojson")
test_df6 <- create_df(geojson_path = "test_gps/20200914101156.geojson")
test_df_list <- list(test_df1, test_df2, test_df3, test_df4, test_df5, test_df6)
test_df_list <- lapply(test_df_list,
function(x) {
names(x) <- c("latitude", "longitude", "Time", "time_long", "time_difference", "speed")
return(x)
})
# Merge lists into one mega dataset
test_df <- do.call("rbind", test_df_list)
# Add combined coordinates to merged sets
test_df$coordinates <- paste(test_df$latitude, ",", test_df$longitude)
# Used this to get the most frequent coordinate, taking latitude and longitude values of each recording as a whole.
# Creating a mode function for coordinates in some sense.
all_coordinates <- table(test_df$coordinates)
test_df$Time <- format(test_df$Time, "%D-%H")
spat_test_df <- SpatialPointsDataFrame(coords = test_df[, c("longitude", "latitude")],
data = test_df[, c("longitude", "latitude", "Time")],
proj4string = CRS("+proj=longlat +datum=WGS84"))
test_df <- test_df[, c("latitude", "longitude", "time_long", "time_difference", "speed")]
test_df
Now, let’s use our function to get our predictions.
assasination_algorithm(test_df)
[1] "Bomb Prediction 1"
[1] "Latitude: 46.8631139796899"
[1] "Longitude: -113.984692535082"
[1] "Time: 1598609011803.84"
[1] "-----------------------------------------------------------------------------------------------------------"
[1] "Bomb Prediction 2"
[1] "Latitude: 46.863371206494"
[1] "Longitude: -113.984826323491"
[1] "Time: 1598607944683.99"
[1] "-----------------------------------------------------------------------------------------------------------"